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1. Introduction 

This book chapter is an eccentric view of the present state of time-dependent 
density functional theory (TDDFT). It is not intended as a comprehensive 
overview of the field, but merely raises some issues that face the field at 
the present time, and we hope it makes enjoyable reading. The opening 
question of Sec. 2.1 is particularly eccentric, mirroring the style of our old 
friend, Bob Parr. 

A time-dependent A^-electron system satisfies the time-dependent Schrod- 
inger equation : 

H9(T 1 ...T N t)=i9(T 1 ...T N t), (1) 

where we have (for simplicity) ignored spin indices, and used atomic units 
(e 2 = h = m = 1), and introduced a dot for time-derivatives. Here the 
Hamiltonian consists of three contributions 

H = f + V ee +V ext , (2) 

the kinetic energy, the Coulomb repulsion, and the external potential, due 
to the nuclei and any external fields. Note Eq. (1) is first-order in time, and 
solutions depend on the initial wavefunction, ^(0). 

Rigorous modern TDDFT begins with the Runge-Gross (RG) theorem 1 , 
although the first modern TDDFT calculations were done by Ando 2 ' 3 
for semiconductor surfaces and by Zangwill and Soven 4 ' 5 for atoms. The 
RG theorem generalizes the Hohenberg-Kohn theorem 6 to time-dependent 
external potentials, and states that, for a given initial state, there is a unique 
mapping between the evolving density and the time-dependent potential. 
We can then consider a system of non-interacting electrons in a Slater 
determinant of orbitals, satisfying: 

| - \ V 2 + v s (rt) J 4>i (rt) = ifc (rt) , (3) 

and beginning in Slater determinant $(0). Their time-dependent density is 



N 

n(rt) = J2\Mrt)\ 2 

i=l 



(4) 
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We can require this density to match that of an interacting electronic 
system, by the RG theorem, and define the time-dependent exchange- 
correlation potential: 

v XG [n; *(0), *(0)](rt) = v s [n; *(0)](rt) - v ext [n; ¥(0)](rt) - v H [n](vt), (5) 

where the Hartree potential is 

«e[ n ](rt) = | d 3 r'^|- (6) 

Note that the exchange-correlation potential is a functional of the initial 
states of both the interacting system and the Kohn-Sham reference system. 
It would be very nice to derive Eq. (3) as a stationary point of an action, 
but this is more subtle than first appears, and is the subject of Sec. 2.2. As 
in the ground-state case, for many electrons, Eq. (3) is much faster to solve 
than Eq. (1). It explicitly yields the time-dependent density n(rt), and, in 
principle, yields everything else we might want to know about the interact- 
ing problem. However, only the density itself is guaranteed to be the same 
in both systems, and even the current, whose gradient is determined by the 
time-derivative of the density via continuity, could differ in the two sys- 
tems, as discussed in Sec. 2.3. Another example is the probability of double 
ionization of an atom in an intense laser field, where the KS expectation 
value differs greatly from the exact value, as shown in Sec. 2.6. Moreover, 
in general, properties of the interacting system are functionals of both the 
density and the initial state, as discussed in Sec. 2.8. 

In practice, we must make approximations for the exchange-correlation 
potential. The most popular in use today is the adiabatic local density 
approximation (ALDA), which employs the potential for a uniform gas of 
electrons of density n(rt): 

«£c LDA (rt)=«"f( n (rt)), (7) 

as was used in the historic calculation of Ando 2 ' 3 . This approximation 
should work well for a system beginning in its ground state which varies 
slowly in time and space, but most of the systems that it is used for do not 
fit this description. The time-dependent Kohn-Sham equations, Eq. (3), 
with the exchange-correlation potential approximated by ALDA, Eq. (7), 
have become increasingly popular for calculations of atoms and molecules 
in strong laser fields, some of which are discussed in Sec. 2.6. How accurate 
such ALDA calculations are has rarely been investigated, but is studied 
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in Sec. 2.7. Note that this approximation is utterly forgetful. The ALDA 
exchange-correlation potential cares nothing for the past history of the den- 
sity or the initial state, and is completely determined by the density at a 
given instant in time. The consequences of this are discussed in Sec. 2.9. 

Most of the applications of TDDFT at present are for the special case 
of an infinitely weak perturbing field applied to a system in its ground 
state. Analysis of this case leads to predictions for electronic transition 
frequencies, oscillator strengths, polarizabilities etc., i.e., the full optical 
response. We define the susceptibility x[no](rr'w) as the response of the 
ground state of the interacting system to a small change in the external 
potential: 



Jn(rcj) = / d 3 r'x[n ](rr'uj)Sv ex t(r'uj) (8) 
dVx s [n ](rr' W )<Mr' W ). (9) 



The second line follows since the density change is the same for the interact- 
ing system and the non-interacting system; x s [no](rr'w) is the Kohn-Sham 
susceptibility. Applying Eq. (5) to a slightly perturbed system, we find 

Sv s (rco) = foext(rw) + j d 3 r' ( j^j^Tj + /xcM(rr'w)^ Sn(r'uj) (10) 

where / X c[^o]( rr 'w) is known as the exchange-correlation kernel. In the 
time-domain, 

f xc [n ](TT l ,t-t')=6v xc [n ](Tt)/6n(T l t l ). (11) 
These equations lead to the Dyson-type response equation 7 

x(rr'cj) =Xs(rr'w) + J d 3 n J d 3 r 2 x s (rriw)/ HXC (rir 2 w)x(r 2 r'w), (12) 

where / HXC = Vl r — r '| + fxc, and all objects are functional of the ground- 
state density. 

Poles of the susceptibility occur at transition frequencies, while oscilla- 
tor strengths are related to pole strengths. Thus solution of Eq. (12) yields 
all optical response information. Casida 8 showed how to solve these equa- 
tions in a finite basis, analogously to the solution of time-dependent Hartree 
equations. Note that explicit calculation requires two separate approxima- 
tions. The first is for the ground-state Kohn Sham system, which yields the 
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ground-state potential, from which is constructed the Kohn-Sham suscep- 
tibility, Xs- The second is an approximation for the exchange-correlation 
kernel, which in principle is frequency-dependent. Sometimes, for v XG a 
standard ground-state functional, such as a generalized gradient approxi- 
mation or a hybrid with exact exchange, is employed, and / xc is approxi- 
mated by the ALDA, i.e., by inserting Eq. (7) in Eq. (11). Other times, both 
are approximated by the same ground-state functional employing, e.g., the 
LDA for v xc and the ALDA for / xc . In this way, the condition 

S 2 E 

lim/ xc (r,r', W )| no = Mr) ^ c (r/) | no (13) 

is satisfied. 

Applications and development of TDDFT is a rapidly expanding field, 
and we next review recent developments. We refer the reader to an earlier 
review 9 and references therein for earlier applications of TDDFT. 

In quantum chemistry, the largest use of TDDFT has been to extract 
the excitation energies and optical response of molecules, using the linear 
response formalism outlined above. Section 2.4 discusses some of the errors 
inherent in such calculations, due to limitations of our present functionals. 
On the other hand, Sec. 2.5 explains how TDDFT can be quite successful 
despite these limitations. 

Most quantum chemical codes, such as Gaussian 10 and ADF 11 ' 12 per- 
form TDDFT response calculations, allowing experimentalists to imme- 
diately compare with theory 13 ' 14 . Surveys of mean polarizabilities of or- 
ganic molecules have been performed 15 . Many molecular calculations, such 
as excitations in small organic molecules 16 , in tertiophene 17 , in transi- 
tion metal molecules 18 and complexes 19 ' 20 ' 21 ' 22 , using non-empirical hybrid 
functionals 23 , have been tested. Closed-shell polycyclic aromatic hydrocar- 
bon cations 24 have been studied, as have been open-shell molecules 25 , lin- 
ear polyene oligimers 26 , radical cations with N-N bonds 27 , and s-tetrazine 
in both gas-phase and solvated 28 . The optical response of organic dyes 29 
has been calculated, and explanations given of the color of 1,2-dithiins 30 , 
and of the spectra of sulfmes 31 , cromone 32 , and pyrazine 33 . Even aro- 
matic radical cations 34 have been studied; also conjugated molecules 35 
and of course fullerenes 36 . Weakly hydrogen-bonded species have also been 
studied 37 . Tests of various functionals for excitation energies on training 
sets of molecules are ongoing 38 , just as in the ground state. One consider- 
able benefit is to combine spectral information with structural and thermo- 
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chemical data as a more thorough test of DFT calculations, as in the case 
of trans-stilbene 39 . 

Higher order response has also been calculated. TDDFT is now used for 
hyperpolarizabilities 40 , such as in chiral molecules 41 . Both linear and mag- 
netic circular dichroism have been calculated 42 ' 43 ' 44 . The Rydberg states of 
propyne have also been calculated in conjunction with resonant-energy mul- 
tiphoton ionization experiments 45 . Recently TDDFT calculations showed 
that a candidate for strong non-linear optical response was less promising 
than experiment suggested 46 . Raman intensities have also been calculated 47 . 
Not yet understood is the ability of TDDFT to perform well, even for some 
states with significant double-excitation character 48 . 

Calculations are appearing with real biological significance, such as on 
chlorophyll A 49 and free-base porphrin 50 ' 51 . A very important recent devel- 
opment is the attempt to study charge transfer in biological molecules 52 ' 53 . 
Fluorescence of 2-aminopurine has been calculated 54 . Algorithmic develop- 
ments are also occurring. Under many circumstances, the Tamm-Dancoff 
approximation is sufficient 55 for excitation energies. Fast algorithms for 
solving TDHF equations can be immediately applied to adiabatic TDDFT 
calculations 56 ' 57 . Geometric derivatives for excited-states have been coded 58 . 

Much of the earliest work on TDDFT response was performed on metal- 
lic clusters 59 . Clusters of ZnS have also been studied 60 . 

In atomic and molecular physics, TDDFT is being applied to problems of 
(mostly) atoms in intense laser fields, including stabilization phenomena 61 , 
and nonsequential multiple ionization 62 ' 63 . Even the more demanding time- 
dependent optimized effective potential (TDOEP) 64 has been coded and 
applied to high harmonic generation 65 ' 66 . Some of these calculations ap- 
ply the basic formalism within Floquet theory 67 ' 68 . Also, the original pho- 
toionization problem has been recently been revisited using more accurate 
ground-state potentials 69 , while calculations for molecules have also been 
done 70 . Clusters can now also be handled 71 . Energetic collisions between 
atoms and ions are also being tackled 72 ' 73 ' 74 . 

For extremely large calculations, involving thousands of electrons, even 
the TDKS equations are too expensive to solve, and time-dependent Thomas- 
Fermi approaches are used 75 ' 76 . 

In mesoscopic physics, TDDFT has been used for quite a while to study 
the optical response of quantum dots 77 . 

The applications mentioned above have all been to finite systems. The 
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extension to extended systems will be a major effort of itself, for the simple 
reason that our present day functionals, such as LDA and GGA, do not 
provide useful approximations to the exchange-correlation kernel / xc in 
this case. To see this, simply consider the Fourier transform of / HX c- The 
direct Coulomb interaction behaves as 47r/g 2 , where q is the wavevector 
corresponding to r — r'. The optical response is the long- wavelength limit, 
dominated by q — > 0. The LDA kernel is local in space, so its Fourier 
transform is a constant, becoming negligible as q — > 0. Similar reasoning 
applies to GGAs. The ultra non-locality in space of the exchange-correlation 
kernel in extended systems has recently been emphasized 78 . 

There have been a few attempts to perform calculations on extended sys- 
tems, such as in polymers 79 ' 80 ' 81 , and the optical response of solids 82 ' 83 ' 84 , 
but these difficulties have not been fully understood or overcome. 

On the other hand, other interesting properties of solids have been 
calculated. For many years TDDFT has been used to calculate the di- 
electric response of metals 85 , in attempts to disentangle band-structure 
effects from correlation, especially in the dispersion of the bulk plasmon 86 . 
Recently, spin-response of magnetic metals 87 ' 88 has been calculated using 
TDDFT. Also electronic damping at surfaces has recently been calculated 
using TDDFT 89 ' 90 . 

At a conceptual level, the links between TDDFT and traditional many- 
body approaches for extended systems, such as GW, are only now being 
explored 91 , as are extensions of ground-state theorems 92 , and the descrip- 
tion of potentials in terms of wavefunction quantities 93 as opposed to func- 
tional derivatives. 

2. Nine questions and some answers 

2.1. Parr question: Do the density and potential determine 
the energy? 

As this volume is dedicated to Bob Parr on his eightieth birthday, we open 
with a question he originally raised for ground-state density functional the- 
ory. The question is: if someone gives you both the external potential for 
a system and it's exact density, can you recover the exact ground-state en- 
ergy, without solving a many-body problem? Note that this differs from the 
usual practical question of DFT, in which you are given only the potential, 
and must find the corresponding density and energy. 
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For two (spin-unpolarized) electrons, the answer is yes. For example, 
one could deduce the ionization energy I from the tail of the density and 
also solve the one-electron problem for the external potential, yielding E\ . 
Then the two-electron ground-state energy is just E\ — I. 

More generally, using a technique invented by Bob and others 94 > 95 , one 
can, for any given density, construct the corresponding Kohn-Sham poten- 
tial, v s (r), i.e., find that single-particle potential for which that density is 
a ground-state density (if it exists). Then one can deduce the exchange- 
correlation potential v xc (r) by inverting its definition in Eq. (5). But note 
that this is not yet enough to determine the exchange-correlation energy, 
just the potential. 

To extract the ground-state energy, we write 



where T s is the Kohn-Sham kinetic energy, U is the Hartree energy, and 
E xc is the exchange-correlation energy. All these pieces can be extracted 
from what we now have, except the last. To go a step further, one can use 
the virial theorem applied to the exchange-correlation potential to find 



where T c is the kinetic contribution to the correlation energy. Thus the 
virial of the exchange-correlation potential yields the sum of energies above, 
but not E xc alone, which would finish our problem. 

Several years ago, one of the authors and Bob (with collaborators) 96 ' 97 
published back-to-back articles in Phys. Rev. A on this point, noting that 
the remaining piece of correlation energy, e.g., E xc — T c , is more amenable 
to approximation by standard density functional methods, so that use of 
the virial above on the exact density reduces errors significantly. 

But, returning to the logical question, note that for high-density sys- 
tems, correlation becomes negligible relative to exchange, so that the answer 
is once again yes. Interestingly, for low-density systems, in which correla- 
tion is dominated by potential contributions, so that T c << \E C \, one can 
use the virial once again to extract E xc , and get the ground-state energy. 

What has this to do with TDDFT? Suppose a system begins in its 
ground-state, and then is disturbed by a time-dependent external potential. 




(14) 




(15) 
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The Heisenberg equation of motion for the Hamiltonian is very simple 92 

E = (dH/dt) = J d 3 r n(rt)v ext (rt) (16) 

i.e., the time-evolution of the energy is given entirely by the time-evolution 
of the one-body perturbation. Integrating with respect to time shows that 
knowledge of v ext (rt) and n(rt) is sufficient to determine the entire evolution 
of the energy, once the initial value is known, i.e., 

E(t) = E(0) + I dt' E(t). (17) 
Jo 

So, rather amusingly, the Parr question is trivially answered yes in TDDFT, 
except for the initial ground-state value. 

So, the Parr question has a positive answer for the ground state of two 
electrons, in the high-density limit, in the low-density limit, and (up to a 
constant) for all time-dependent problems. Is it true in general? 



2.2. What is the simplest definition of the action? 

Time-independent quantum mechanics is armed with a variational principle 
which is very useful for ground-state DFT. By approximating the exchange- 
correlation energy functional, we obtain approximations for its functional 
derivative, the exchange-correlation potential w X c[^]( r ) to be used in Kohn- 
Sham calculations. 

What is the analog in the time-dependent case? In time-dependent quan- 
tum mechanics, the role of the energy is taken by the quantum mechanical 
action: 

A[9] = dt(9(t)\idt - H(t)|*(t)>. (18) 

Stationary points of A, with initial wavefunction *&(to), yield solutions to 
the time-dependent Schrodinger equation. Translating this into TDDFT is 
however much more subtle than in the ground-state case. The simplest, and 
perhaps most natural step would be to define 

A[n; * ] = / 1 dt(9[n; *o](*)|»St " H 0o (t)|¥[n; * ](*)>, (19) 

where *o] is a functional of the density and the initial state *o 1 - 
The potential vq in the Hamiltonian H Vo is given and fixed. The Euler 
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equation 6A/6n(rt) = would then supposedly yield the correct density. 
One can extract out a universal part of the action 98 ' 9 , and a part that 
explicitly depends on the external potential. One can define a similar object 
for the Kohn-Sham action 98 ' 9 and, comparing the two, extract an exchange- 
correlation action .Axcfa; ^o, $o]- 

There are, however, a number of problems with this definition of the 
action. First of all the density n and the initial state determine *f>o](t)) 
only up to an arbitrary phase factor. This means that the wavefunction 
|$(*)) = exp(ia(t))\^[n;^ ](t)) with a(t ) = but otherwise arbitrary, 
also gives the density n and the initial state \^(to)) = |*o)- This makes the 
action ill-defined since 



<*[n; *o](*)|»St - H(*)|*[n; *o](*)> = - #(*)!*(*)> + d t a(t). (20) 



It is easily seen that the effect of a(t) is to add a purely time-dependent 
shift d t a(t) to the potential that |\P(£)) evolves in. To make the action well- 
defined we have to specify a(t), which means that we have to fix a gauge for 
the potential. One choice would be that we choose as an argument of the 
action functional the wavefunction *o](*)) that evolves in a potential 
v(rt) that vanishes at infinity. This defines the potential uniquely in terms of 
the density and the initial state. If we make this choice then we can express 
the action explicitly in terms of the density n, the potential v[n;^o](rt) 
and the potential vo(rt) of the Hamiltonian that defines the action: 



where we used that *o](*)) satisfies the time-dependent Schrodinger 

equation (TDSE) with potential v[n,^o](rt). If we calculate the functional 
derivative of this action with respect to n we obtain 



A[n; *o] = / 1 dt(9[n; *o](*)|»St " H„ (t)|*[n; * ](*)> 



dt{9[n; *o](*)|»St " H v (t) + V[n; * ](*) - H(i)|*[n; * ](*)> 



to 




(21) 



Sn(rt) 



SA 



= v 




Svjr't 1 ) 
Sn(rt) ' 



(22) 



This is however not what one at first sight would have expected. One might 
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have expected that 



SA 



v[n,* ](rt) - v (rt) 



(23) 



Sn(rt) 



which makes the action stationary for v[n, ^o](rt) = vo(rt); then the wave- 
function *o](*)) which makes the action stationary is the one that 
satisfies the TDSE with Hamiltonian H Vo . This is, however, not the case. 
So what has happened? To understand this point we must go back to the 
original action Eq.(18) as a functional of the wavefunction, rather than the 
density, and see under which conditions we can derive the TDSE from the 
action. In fact, only under certain types of variations do the stationary 
points of the action yield the TDSE (id t — H)\W) = 0. We refer the reader 
to the recent review 99 for details, and here just state the results. Suppose 
the action A is stationary for variations around a certain ty. Then 
satisfies the TDSE if the variations are such that 

(1) S^(to) = S^(ti) = and the real and imaginary part can be varied 
independently, 

(2) or, alternatively, both = <5$ and 6^2 = are allowed variations 
for any <5$ 100 . 

Let us now go back to our TDDFT action of Eq.(19). This action is obvi- 
ously defined on a restricted set of wavefunctions, namely all wavefunctions 
that can be parameterized by densities on the basis of the Runge-Gross 
theorem together with the gauge condition v — > for |r| — > oo on the corre- 
sponding potential. We will call this set of wavefunctions V. To see if we can 
derive the TDSE from the restricted TDDFT action we must check that 
within the restricted set of wavefunctions V we can make the variations 
mentioned in points (1) and (2) above. We will see that this is not possible 
which is then consistent with Eq.(22). The variations 6V within the set V 
must always be generated by potential variations, i.e. if some ^ in this set 
satisfies a TDSE with Hamiltonian H then \l> + ^\l> satisfies 



(idt-H v -6V)\V + 6V) =0 



(24) 



for some potential variation SV . If we collect the first order terms we see 
that 6V must satisfy 



(idt - H v )\89) = 8V\9) 



(25) 
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with the boundary condition S^(to) = since we evolve all wavefunctions 
in the set V from a fixed initial state iPo- Now it is clear that the real and 
imaginary part of are not independent, they are both determined by 
the potential SV . Moreover, we see that Eq.(25) is first order in time and 
therefore S^f(t) is completely determined by the initial condition S^(to) = 
0. We are therefore not allowed to put a second constraint S^(t\) = on 
the variation S^. Therefore we can not do the variations mentioned in point 
(1) above and are therefore not able to derive the TDSE from the TDDFT 
action in this way. What about the variations mentioned in point (2)? 
We see immediately that if S^i = <5$ is a variation that satisfies Eq.(25) 
then the variation 6^2 = is produced by the potential iSV, which 
is imaginary and therefore not allowed as a potential. We can therefore 
not make the variations mentioned in point (2) and we see that we can 
not derive the TDSE from the TDDFT action in this way either. We thus 
conclude that the action can not be used as a basis of time-dependent 
density functional theory. 

The obvious question is then, can we define some other action functional as 
the basis of a time-dependent density functional theory? For instance, can 
we find some action A that satisfies Eq.(23)? If we can find such an action 
then we can construct the Legendre transform 

A[v] = -A[n] + J d 3 rdtn(rt)(v(rt) - v (rt)), (26) 

where the density n on the right hand side must now be regarded as a 
functional of the potential v. The functional A[v] satisfies 

S A f SASn f. .Sn 

— = - ——+ (v-v )—+n = n, (27) 
ov J on ov J ov 

where for convenience we used a shortened notation and left out the argu- 
ments. We see that if we can find a functional A[v] of the external potential 
such that 

SA 

n(rt) (28) 



Sv(rt) 

then the functional A[n] with property Eq.(23) can be constructed from 
the inverse Legendre transform: 

A[n] = -A[v] + J d 3 rdtn(rt)(v(rt) - v (rt)), (29) 
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where v must now be regarded as a functional of n. The question therefore 
is: Is there a functional A such that Eq.(28) is satisfied? We therefore ask 
the question whether or not the fundamental variable of TDDFT, namely 
the time-dependent density n, can be obtained as the functional derivative 
of some functional A of the external potential. The answer is no if A is 
twice differentiable because then we could differentiate Eq.(28) again and 
obtain 

S A = Wrt. (30 ) 

6v(rt)6v(r't') 6v(r't')' V ' 

The left hand side of this equation is symmetric in the space-time arguments 
since we assumed that A was twice differentiable, whereas the right hand 
side of this equation is the density response function which has a causal 
structure, i.e. it is zero for t > t' . Therefore the causality and symmetry 
requirements contradict each other. We conclude that there is no differ- 
entiable functional of the external field with the property Eq.(28). Con- 
sequently there is no functional of the density with the property Eq.(23). 
We therefore conclude that the external potential v[n, ^o] of a many-body 
system with initial state and density n can not be obtained as the deriva- 
tive of a density functional. The same is of course true for a noninteracting 
system and in particular the Kohn-Sham system. We must conclude that 
the time-dependent Kohn-Sham potential is not a density derivative. To 
be more precise, there is no functional of the w-representable density n(rt) 
that has the Kohn-Sham potential as its derivative. 

Is there some action functional defined on a larger class of densities that is 
capable of providing a derivation of the time-dependent Kohn-Sham equa- 
tions? The answer to this question is yes. It turns out that that one can 
define a functional on a set of so-called time contour densities from which we 
can derive the Kohn-Sham equations. The corresponding action functional 
is called the Keldysh action. For further details on this functional we refer 
to the literature 101 >". Let us finally answer the topical question asked in 
the title of this section. What is the simplest definition of the action? At the 
moment this is certainly the Keldysh action, as this is currently the only 
action functional that leads to a derivation of the Kohn-Sham equations 
that is free of paradoxes. 
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2.3. Is the Kohn-Sham current equal to the true current? 

In TDDFT a noninteracting Kohn-Sham system is introduced with the 
same density as the true interacting system. Since the density and the 
current are closely related by the continuity equation we may wonder what 
the relation is between the Kohn-Sham current and the true current. In 
order to avoid confusion we stress that we will always be dealing with 
systems in time-dependent external fields that can always be transformed 
to a pure scalar potential by a gauge transformation. We therefore exclude 
magnetic fields. We start by describing some properties of the current. The 
current is defined as 

j(rt) = ^(V-V , )7(r,r , ,t)| r=r . (31) 

where 7 is the one-particle density matrix of the system. The expectation 
value of the momentum P can be directly calculated from the current as 
follows 

P(t) = J d 3 rj(rt). (32) 

Finally, if we calculate the commutator of the density operator with the 
Hamiltonian we obtain the continuity equation 

Stn(rt) = -V-j(rt). (33) 

Let us now turn to the Kohn-Sham system which has a current 

1 N 

j*( r *) = ^ E (v*(rt)V^(rt) - (V<pt(Tt))<p k (Tt)) (34) 

and where ipk are the Kohn-Sham orbitals. The Kohn-Sham current is not 
necessarily equal to the true current. We therefore define the exchange- 
correlation part j xc of the current by 

j(rt)=j s (rt)+j ac (rt). (35) 

Let us now investigate what the equations we defined at the beginning of 
this section can tell us about j xc . We start by considering the continuity 
equation Eq.(33). Since the true system and the Kohn-Sham system by 
definition have the same density we obtain 



V ■ j ac (rt) = V ■ j(rt) - V • j,(rt) = 0. 



(36) 
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We see that j xc is a divergenceless or tranverse vector field which can there- 
fore be written in the form j xc = V x C for some vector field C. The con- 
tinuity equation has another consequence. If we consider finite system for 
which densities and currents vanish at infinity then using Eq.(33) we can 
write the expectation value of the momentum as 



From this equation we see that the Kohn-Sham system and the true system 
have the same momentum and so 



The continuity equation has therefore told us that j xc is a transverse vector 
field whose spatial average is zero. So far we have discussed some general 
relations which must be satisfied by j xc . We now turn to some more specific 
cases. For one-dimensional systems in which there are no transverse vector 
fields j xc vanishes. For such systems the current is uniquely determined by 
the density from 



if the current vanishes at infinity. This equation immediately implies that 
for one-dimensional systems j = j s and j xc = 0. It was conjectured 98 that 
this holds true in general. The argument was based on the fact that systems 
with two different time-dependent scalar potentials v ^ v' + C(t) differing 
by more than a purely time-dependent function, yield two different currents 
j ^ j'. This follows immediately from the proof of the Runge-Gross theo- 
rem. Therefore the external potential is a well-defined functional on 
the set of w-representable currents. However, it was implicitly assumed that 
every current j produced by a scalar potential v in an interacting system 
can also be produced by a scalar potential v s in a noninteracting system, 
i.e. the proof was based on an unproven noninteracting w-representability 
assumption. One may object that such an assumption is also made for the 
density in constructing the Kohn-Sham system in the first place. However, 
the density n and the scalar potential v are conjugate variables in the sense 
that the contribution of the external potential to the total energy is of the 
form of an integral over n times v. It is exactly this property that is used in 
the proof of the Hohenberg-Kohn theorem. The conjugate variable of the 




(37) 




(38) 




— CO 



(39) 
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current, on the other hand, is the vector potential. On the basis of this one 
would expect that to reproduce an interacting current in an noninteract- 
ing system one needs an exchange-correlation vector potential A xc . Such 
a vector potential is indeed introduced in time-dependent current-density 
functional theory. On the basis of these arguments it seems unlikely that 
the noninteracting w-representability assumption for the current can be jus- 
tified. 

Another argument which suggests that j xc is nonzero in general is provided 
by the following example. Consider a Kohn-Sham system with two particles 
in a singlet state, doubly occupying one spatial orbital: 

id t ip(rt) = ("^V 2 + v,(rt)j <p{rt) 

n(rt) = 2\ip(rt)\ 2 . 
The Kohn-Sham orbital can then always be written in the form 



V>(rt) = y^e*< rt >. (40) 

With this expression we obtain the following equation for the Kohn-Sham 
current 



j,(rt) = n(rt)V0(rt) (41) 



and we see that 



V x = 0. (42) 

n(rt) v ' 

If we assume that j = j s then this implies that for any interacting two- 
electron system 

Vx^=0. (43) 
n(rt) 

This seems an unlikely property for an arbitrary two-electron system with 
rotating currents, such as an Helium atom in an intense laser pulse of circu- 
larly polarized light. If the property Eq.(43) would be a general feature of 
two-electron systems then it must follow from some special property of the 
two-particle Hamiltonian. The question would of course be settled with one 
counterexample for which the vorticity of Eq.(43) does not vanish. How- 
ever, it is not easy to give a simple example. The solvable model systems 
are separable and have the special feature that they decouple relative from 
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the center-of-mass motion which leads to currents for which Eq.(43) seems 
to be true. For instance the separable harmonic atom has a rigid mode 
with current j(ri) = n(ri)dR/di where R(i) is the expectation value of 
the center-of-mass of the two-electron wavefunction. We therefore have to 
wait for numerical counterexamples of nonseparable systems, which must 
at least be two-dimensional. If such an example of nonvanishing vorticity 
can be found it would imply that, generally, the true current density of 
an interacting two-electron system cannot be reproduced by the current 
density of a spin restricted Kohn-Sham system. However, if we relaxed the 
restriction that the Kohn-Sham system is in a doubly-occupied spatial or- 
bital, and instead considered singlet states with two distinct orbitals, then 
Eq. (42) does not necessarily hold. We may ask, for a given time-dependent 
density and current, can we find a Kohn-Sham system which reproduces 
both the interacting density and interacting current? If so, is this Kohn- 
Sham system unique? In one-dimension, the answer to the latter question 
is, in many cases, no, but what about in more than one dimension? 

2.4. What errors does my TDDFT calculation of electronic 
transitions make? 

As mentioned in the introduction, there are two approximations in any prac- 
tical TDDFT calculation of transitions: the approximation for the ground- 
state potential and that for the XC kernel. 

To study the first source of error, we begin with atoms, and then move 
on to molecules and solids. Interestingly, most of our presently used ground- 
state XC energy functional have potentials that do not resemble the exact 
XC potential very closely, especially for the He atom. In Fig. 1, we plot the 
XC potential for LDA, the PW91 GGA, and the exact potential. In the 
asymptotic region, the exact exchange potential decays as — 1/r, whereas 
the LDA and GGA decay exponentially. (How such functionals still yield 
good ground-state energies is an interesting question in itself 102 ). This is 
not a difficulty for the ground-state theory, as the potential in the region 
where bound orbitals (in this case, just one) live is well-approximated. But 
excited states, especially Rydberg states, are badly described. In fact, most 
of them are unbound. 

How can we avoid this problem? The least expensive approach is to fix 
up the potential by hand, by adding the known —1/r tail in the asymp- 
totic region. Several prescriptions for doing this have been suggested in 
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Fig. 1. The exchange-correlation potential for helium: exact (solid), exact exchange- 
only approximation(dotted), local density approximation (dashed) and the PW91 gen- 
eralized gradient approximation (dashed-dotted). 

the literature 103 ' 104 ' 105 . A statistical average of orbitals also produces an 
excellent potential in this regard 106 ' 107 . Alternatively, the technology for 
including exact exchange in DFT calculations has been greatly developed 
over the last decade due to the advent of the KLI approximation, so that 
it has become possible for molecules. The exact exchange potential cannot 
be distinguished from the exchange-correlation potential in Fig. 1 and the 
differences in unoccupied orbital energies and matrix elements over these 
orbitals are tiny 108 ' 109 ' 110 . These problems have recently been reviewed by 
Tozer and Handy 111 . There is a noticeable error still in the occupied orbital 
energy, since the correlation potential is more significant near the nucleus. 
Hybrid functionals, which mix in only a fraction of exact exchange with 
GGA, only partially cure the problem 112 . 

These problems are most dramatic in the long-range decay of the po- 
tential, which dominates all unoccupied states for the He atom. As we 
consider larger systems, the problem gets less. Even for the Be atom, the 
first transition from 2s to 2p is reasonable in LDA or GGA, but higher lev- 
els are bad. For many molecules and reactions of photochemical interest, it 
is only excitations to the first few low-lying states that are important, and 
so these asymptotic difficulties are irrelevant. In the limit of bulk solids, 
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these difficulties do not exist. 

Of course, if the underlying ground-state problem is not weakly cor- 
related, then errors in the ground-state approximation become large, and 
difficulties arise for excitations. The H 2 molecule, stretched beyond the 
Coulson-Fisher point, so that an LDA or GGA ground-state calculation 
spontaneously breaks symmetry, is a very demanding case for TDDFT, 
being based on the KS determinant, and has been recently studied from 
several different angles 113 ' 114 ' 115 . 

Assuming the ground-state potential is accurate, how about approxi- 
mations to /xc? In almost all applications at present, / xc is approximated 
adiabatically, by the second derivative of a ground-state XC energy func- 
tional. In recent work 116 , some of us investigated the effect of different 
approximations. Using the single-pole approximation (discussed in the next 
section), we analyzed the results, finding that they could be understood on 
the basis of trends already known for ground-state functional . We found 
that the mean of the singlet and triplet levels of the He atom differed from 
the exact KS transition by only an expectation value of the parallel-spin 
correlation contribution in / xc , typically a very small part of E xc . This 
gives a functional explanation of why KS values are often good zero-order 
approximations (see the next section for more detail). We also found that 
the splitting depended on a cancellation between exact exchange and an- 
tiparallel correlation. Thus a hybrid of exact exchange for the mean and 
ALDA for splitting led to very good agreement with experiment for the 
spectrum of the He atom. This hybrid was designed simply to illustrate 
how insight into functionals could be used to improve accuracy in TDDFT 
calculations. Its construction depended on the system being weakly corre- 
lated and having few electrons. It worked for Be, but less well, and has been 
shown to yield no improvement for stretched HHe + , or for the dispersion 
of the plasmon in the uniform electron gas. 

Very little exact information is known about the kernel. Even for exact 
exchange, relatively little practical is known 117 , beyond the two-electron 
unpolarized case. The frequency-dependence of / xc for the uniform gas has 
been under constant study 118 ' 119 ' 120 ' 121 ' 122 - 123 . 

As for the question of which approximation - the one for v xc or the one 
for /xc - has a stronger influence on the calculated spectra, generally the 
effect of v xc is much stronger for higher-lying excitations. A typical picture 
is shown in Fig. 2, where we compare the errors of the singlet excitation 
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Fig. 2. Errors of singlet excitation energies from the ground state of Be, calculated 
using the exact exchange-correlation potential obtained from an accurate wavefunction 
calculation, the OEP-SIC approximation and the X-only KLI approximation and with 
different approximations for the exchange-correlation kernel. The errors are given in 
mHartrees. To guide the eye, the errors of the discrete excitation energies were connected 
with lines. Taken from Ref. 116 



energies of the Be atom resulting from various approximations for v xc and 
/xc- For low-lying excitations, and especially for larger systems, often it 
is /xc which has the larger effect, as is seen in the 2s2p transition in the 
figure. These lower lying excitations are often the ones of photo-chemical 
interest. 



2.5. When are Kohn-Sham transitions good 
approximations ? 

Early on, Casida 8 showed how to recast Eq. (12) into common quantum 
chemical notation 



KF = n 2 F. 



(44) 
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The eigenvalues are the squares of the physical excitation energies ft, and 
the eigenvectors F determine the true oscillator strength of the correspond- 
ing transitions 8 . Alternatively, these equations may be written 108 : 

^2(M qq ,((l)+u q S qq .)P q . =(ip q . (45) 
q' 

Here the matrix elements M qql (ui) = J d 3 r J dV$*(r)/ HXC (r, r', cj)$,<(r') 
are expressed in terms of Kohn-Sham orbitals $ g (r) = ^ a (r)^»,(r). The oj q 
denote Kohn-Sham transitions q = i,a from an occupied KS orbital i to an 
unoccupied KS orbital a. 

To gain more insight into the structure of the solutions to Eqs. (45) and 
(44), we first notice that, if / xc is zero, the matrix is diagonal, and both 
eigenvalues and eigenvectors will equal their KS counterparts. For weakly 
correlated systems, we expect / xc to be, in some sense, small, so that the 
KS values should be good approximations. But, small compared to what, 
i.e., when will the corrections to KS values be accurate? 

To analyze these corrections, we assume for simplicity that / xc is frequen- 
cy-independent, as is the case in most currently used approximations. We 
then use the method of continued fractions (CF) 124 , where solutions of 
linear equations are given in the form of a continued fraction expansion. 
Truncating the continued fraction at a given order corresponds to perform- 
ing perturbation theory in the distance from the diagonal of the matrix. 
This is very different from the traditional Gorling-Levy (GL) perturbation 
theory 125 > 126 5 which is an expansion in powers of the adiabatic coupling 
constant A. The CF expansion does not assume / xc is small, and can be 
valid even when there are large corrections to KS results. Starting from a 
continued fraction expression the conventional perturbation expansion may 
be recovered by a consistent Taylor expansion in the coupling parameter 
up to the desired order. Using such a low order truncated CF expansion we 
find for the excitation energies in our eigenvalue problem 

9 99 ^ {u q + M qq ) - (u) r + M rr ) 
y* M qr M rs M sq | 

The asterisk on the sums indicates that the summation is only performed 
over terms having distinct indices. (At this point we should note that our 
result (46) remains valid even for frequency dependent matrix elements. In 
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this case a truncated version of (46) becomes a nonlinear equation in Q.) 
Looking at the result, we see that the CF method produces a power series 
in a (hopefully) small parameter, which is a matrix element divided by a 
transition frequency difference. Thus, truncation of this series is accurate 
when this ratio is small, i.e, when the shift away from KS transition fre- 
quencies is small relative to the separation. It is an approximation of weakly 
coupled transitions, rather than of weakly correlated systems. 

The zeroth order result was dubbed the single-pole approximation (SPA) 
7 , and was used in early calculations of transition frequencies. It has also 
been used to relate TDDFT results to GL perturbation theory results for 
excitations 91 . But the expression in Eq. (46) tells us more. Clearly, the SPA 
will be a good approximation to the full excitation energies only if higher 
order terms like the sums in (Eq. 46) contribute little. We can expect the 
SPA to be valid if the SPA shift of a transition is small relative to the 
separation between that transition and its neighbors. By expanding / xc 
to second-order in A, it also yields an exact expression for the transition 
frequency to second-order in GL perturbation theory. 

Let us now take a look at the eigenvectors. For this purpose we use the 
second version of the eigenvalue problem Eq. (44). Applying the CF method 
to Casida's formalism, we find the oscillator strength expansion to be 



Here co q and fi q denote KS excitation energies and dipole matrix elements 
respectively. The leading term is the SPA result, and is simply the KS oscil- 
lator strength. Thus, contrary to the transition frequency, there is no cor- 
rection to oscillator strengths in SPA. This can be easily seen from Eq. (46). 
If only diagonal matrix elements are retained, the eigenvectors remain unit 
vectors and don't change. But once again, with correction terms, we un- 
derstand much more. Corrections to KS oscillator strengths remain small 
if transitions are well-separated, even for strongly correlated systems. Also, 
Eq. (47) is the exact GL expansion for oscillator strengths, if / x is used for 
the matrix elements. 

We end this section by studying the effect of approximate ground-state 
potentials on oscillator strengths. As in orbital energies, the potential must 
have a long-range decay to bind Rydberg states. Once that is so, the KS 
oscillator strengths depend on both transition frequencies and dipole matrix 




(47) 
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elements. In numerical studies on the He atom, we find that the dipole 
matrix elements in exact X or LDA-SIC potentials are extremely close to 
the values in the exact XC potential 127 , and that the main source of error in 
such approximate calculations is then in the orbital energies, as discussed 
in the previous section. 

2.6. What happens to atoms in strong laser fields? 

TDDFT really takes off (as do the dynamics) when atoms and molecules are 
subjected to strong laser fields. Many new phenomena such as multiphoton 
ionization (MPI), 128 - 129 - 130 - 131 above-threshold ionization (ATI) 132 . 133 . 134 
or high-harmonic generation (HHG) 135 > 136 are observed when the electric 
field amplitude of the laser is comparable to or even exceeds the static nu- 
clear Coulomb field experienced by the electrons. TDDFT is perhaps the 
only feasible method to calculate the time-dynamics of interacting many- 
body systems in this regime. In spite of the fact that the electron-electron 
interaction is much weaker than the strong external driving field, electron 
correlation effects can be important 137 > 13 M39,i40 F or anv time-dependent 
calculation within TDDFT there is, however, one important point to note: it 
is not sufficient to know a good approximation for the exchange-correlation 
potential (and thus a good approximation for the time evolution of the den- 
sity) . It is in addition necessary to know density-functionals for the observ- 
ables. Usually the quantities of interest in a time-dependent calculation are 
some generalized cross sections like S or T matrices in scattering processes, 
ionization yields, branching ratios for chemical reactions or dissociation 
probabilities, to mention just a few. By virtue of the Runge-Gross theorem 
all of these observables are functionals of the time-dependent density and 
the initial-state. Consequently, a calculation within TDDFT involves two 
steps 

(i) The TDKS equations are solved using some approximate form for the 
exchange-correlation potential. The density is calculated from the time- 
dependent orbitals. 

(ii) The (approximate) time-dependent density from step (i) is then in- 
serted in the functionals for the physical observables of interest. 

In some cases we know the exact functional dependence of such generalized 
cross sections on the time-dependent density, but in most cases we do not. 
One prominent example where the exact form of the functional is known is 
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the case of harmonic spectra for atoms or molecules. Neglecting propaga- 
tion effects in the medium, as well as the fact that the focal volume has an 
intensity profile exposing different atoms to different intensities, these spec- 
tra are given by |d(w)| 2 . Here d(oj) is the Fourier transform of the induced 
time-dependent dipole moment of the system 

d{t) = J d 3 rzn(rt). (48) 

Ionization yields and photoelectron spectra are much harder to express as 
functionals of the density. For the ionization yields of the Helium atom 
some approximate functionals are known 141 . The construction of function- 
als yielding probabilities P +1 and P +2 for singly and doubly ionized Helium 
rests on a geometrical concept. Consider the spatial partitioning of the norm 
of the true time-dependent two-particle wavefunction in the following form 



1 



/ d 3 n [ d 3 r 2 |*( ri ,r 2 ,i)| 2 +2 / d 3 n I d 3 r 2 |*( ri ,r 2 ,i)| 2 

J A J A J A JB 

+ I d 3 n I d 3 r 2 |*( ri ,r 2 ,i)| 2 . (49) 

JB JB 

Here the analyzing volume A describes an appropriately chosen region 
which encloses the nucleus and B is its complement B = R 3 \ A. The 
factor of two in front of the second integral appears because of the anti- 
symmetry of the wavefunction. Due to the probability interpretation of the 
wavefunction the second integral (AB) in Eq. (49) is equal to the probabil- 
ity of finding one electron inside the volume A and simultaneously finding 
a second electron outside the volume A. This can be viewed as single- 
ionization P +1 . Similarly the integral (BB) corresponds to the probability 
of double-ionization P +2 . Introducing the pair-correlation function 

fl(r 1 ,r 2 ,t)= , r(r .\' r2 ; f) . v (50) 

where T(ri,r 2 ,i) is the time-dependent two particle density matrix (in 
the case of Helium simply T(ri,r 2 ,i) = 2|$(ri,r 2 ,i)| 2 ) the expressions 
for single- and double-ionization of Helium can be written as 

P+\t) = f d 3 rn{r,t)- f d 3 n f d 3 r 2 n(ri , t) n(r 2 , t) g[n]{r u r 2 , t) 

J A J A J A 

P+' 2 (t) = l- [ d 3 rn(r,t) + l [ d 3 n [ d 3 r 2 n(rj , t) n(r 2 , t) g[n]{r u r 2 , t). 

J A 1 J A J A 



August 14, 2001 8:1 WSPC/Guidelines n8'13 



Ten topical questions in time- dependent density functional theory 25 

(51) 

Note, since the time propagations are started in the ground state, both ion- 
ization yields P +1 ,P +2 are functionals of the density only. For Helium, the 
simplest approximation for the pair-correlation function is the exchange- 
only expression 

flx[n](r 1 ,r 2 ,t) = i. (52) 

In this case the ionization probabilities P +1 ,P +2 reduce to 
P + \t)=2p{t){l-p{t)) 

P +2 (t) = (l-p(t)) 2 , (53) 

where p(t) = \ f A d 3 rn(r,t). Note, that exactly the same result (53) is 
recovered when a product wavefunction is inserted in Eqs. (51). 

To assess the quality of the approximation (52) involved in the function- 
als P +1 ,P +2 Lappas and van Leeuwen 144 have performed numerically exact 
time propagations for a ID soft-core model of Helium in a laser field 145 . In 
the length gauge the Hamiltonian for their model system reads 

+(xi + x 2 ) E f(t) sin(cji). (54) 

Here f(t) describes the envelope of the laser pulse and V(x) = l/\/x 2 + 1 is 
the soft-core model potential. From the time-evolution of the two-electron 
wavefunction, using the split-operator method 146 , numerically exact ref- 
erence data for the ionization yields were obtained from expressions (51). 
This was then compared to ionization yields from Eq. (53) using the exact 
density calculated from the correlated wavefunction. In this way approxima- 
tions in the first step (i) of the computational procedure were circumvented 
and the accuracy of the functionals for the cross sections can be tested 
directly. 

The results of Lappas and van Leeuwen are shown in Fig. 3 where the 
ion yield for single (triangles) and double (squares) ionization is plotted 
as function of the laser intensity. Although the double-ionization evaluated 
from the approximate functional (53) still does not agree with the exact 
yield, the famous knee structure 147 ' 148 is reproduced. This is in contrast to 
TD x-only calculations where the knee structure cannot be recovered at all. 
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Fig. 3. Single- and double-ionization yields of He from the fully correlated exact model 
(full triangles and squares), and the yields based on Eqs. (53) in the text evaluated with 
the 'exact' electron densities (open triangles and squares). Taken from Ref. 144 



In our context this situation would amount to the evaluation of (53) with 
the approximate TD x-only density. 

In summary we conclude that the two step procedure (i),(ii) is necessary 
for any time-dependent calculation in TDDFT. Both an approximation 
for the exchange-correlation potential and a functional approximation for 
the generalized cross section of interest, have to be known. The relative 
importance of these different types of approximations has to be investigated 
for any particular case of interest. 

In the special case of Helium double-ionization it turned out that the 
functional approximation for the ion yields (53) recovers the well known 
knee structure only when the functional were evaluated with the exact 
densities. Evaluation with approximate TD x-only densities did not repro- 
duce the essential physics 141 ' 142 . Thus more accurate exchange-correlation 
potentials, possibly including memory effects 143 , have to be utilized in this 
case to obtain better approximations to the true time-evolving density. 
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2.7. When does ALDA work beyond the linear response 
regime ? 

Most of the TDDFT calculations use the simple ALDA approximation 
which treats the instantaneous density as if it was a ground-state den- 
sity and applies the local density approximation for the exchange correla- 
tion potential Eq.(7). LDA is ubiquitous and reliable for most ground-state 
systems, but the time-evolving system is certainly not typically a ground- 
state, and memory effects are neglected (see Sec. 2.9). How well does such 
an approximation work? In the ground-state case, conditions that exact 
functionals satisfy are an important guide to understanding why various 
approximations work (or fail) as well as they do, and indeed, to their con- 
struction 149 . Exact conditions that the time-dependent functionals satisfy 
include Newton's Third Law 15 °. 151 . 152 5 the harmonic potential theorem 153 , 
behavior under uniform scaling 92 , a virial theorem 92 , and the memory for- 
mula 154 . For example, the violation of the harmonic potential theorem by 
the Gross-Kohn approximation 123 for the exchange correlation kernel pro- 
vided much of the motivation for a search for other approximations which 
do satisfy this theorem 155 . 

These exact conditions provide only a small number of tools, albeit 
important ones, to take in the knapsack when exploring the vast expanse 
of possible dynamical behavior. Much has yet to be learnt about properties 
of time-dependent functionals in order to obtain accurate approximations. 

A large part of the problem is that until two years ago, there were no 
exact time-dependent Kohn-Sham calculations done; that is, exact calcula- 
tions of an interacting system and of the corresponding Kohn-Sham wave- 
function. For this purpose, time-dependent Hooke's atom, two interacting 
electrons in a harmonic well of time-dependent force constant has been ex- 
ploited in three recent works 156 > 92 > 157 : these are the first numerically exact 
time-dependent Kohn-Sham calculations of any system. 

The two electrons in time-dependent Hooke's atom live in the Hamilto- 
nian 

H = ~\ (V 2 + V 2 ) + \k{t) (r 2 +rl) + (55) 

where the time-dependent force constant is 

k(t) = k — ecos(wi). (56) 
The one-electron version of this is the Mathieu oscillator 158 . The dynam- 
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ics is largely classical, because of the quadratic nature of the potential, and 
classically evolving a bunch of trajectories with the same initial phase-space 
distribution as the initial ground-state, describes the full quantum dynam- 
ics well. For e not too large, time-dependent perturbation theory tells us 
what frequencies appear in the dynamics. In Fig. 4, we plot the root-mean- 
square variance of the distribution as a function of time for one electron, 
calculated quantum mechanically, quasi-classically and in perturbation the- 
ory. The rough trends in the two-electron density are similar to those of the 
one-electron density, and this is also shown here. The difference is due to 
interaction. 

J<r 2 >(t) 




1.6 - 



I I I I I I I I I L 

3.6 I 1 1 1 1 1 1 1 1 r 




1.4 b i i i i i i i l 

5 10 15 20 25 30 35 40 45 
t 



Fig. 4. Rms variance: one electron in the Mathieu oscillator, calculated quantum me- 
chanically (solid line in bottom panel); time-dependent perturbation theory (dashed line 
in bottom panel); quasiclassically (dotted line); two electrons of time-dependent Hooke's 
atom (middle panel); the ground-state corresponding to the external force constant (top 
panel). The parameters in Eq.56 are k = 0.25, ui = 0.75 and e = 0.05. 

The two electron case may be solved exactly numerically. Transforming 
to center-of-mass and relative coordinates renders the Hamiltonian sep- 
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arable, and thanks to spherical symmetry, one needs only to solve two 
uncoupled one-dimensional time-dependent Schrodinger equations numeri- 
cally. The calculation begins in the ground-state. Thus the exact evolving 
wavefunction and density are obtained. 

Since there is just one (complex) occupied orbital in the Kohn-Sham 
calculation, a simple inversion of the Kohn-Sham equation yields the Kohn- 
Sham potential in terms of the evolving density. Once v s (r,t) is known, the 
energy components T s (t) ,U(t), E xc (t) may then be extracted with the help 
of the equations of motion 92 . These energies contain global information 
about the potentials v s (rt) , v xc ( T t) . We shall describe a manifestation of 
this shortly. 

The exact calculation may be compared with that of an exact adiabatic 
calculation 157 . To this end, we define the ground-state components of the 
various energy components as the value of the exact ground-state functional 
evaluated on the instantaneous density. The difference between this and 
the exact energy is termed the "dynamical component", e.g. Et vn {t) = 
E c (t) — Ec S [n(rt)]. The exact ground-state value is obtained by observing 
that, for our choice of time-dependent potential, the static potential which 
has the ground-state density matching the instantaneous density at time 
t is very close to that of a static Hooke's atom of a certain force constant 

k ef f{t). 

Dynamical effects were found to be very large: except when the force 
constant is varied slowly enough that the system remains in the instan- 
taneous ground-state, the functionals behave qualitatively differently than 
the adiabatic approximation. We refer the reader to the paper 157 for many 
interesting results from a variety of runs and here just present a few. 

In Fig. 5 we plot the correlation energy and its first time-derivative for a 
typical run. A feature which would put a ground-stater out the door is that 
E c (t) can become positive. No adiabatic approximation can capture this. In 
the time-dependent case there is of course no variational principle holding 
E G down below zero and we found positivity in all of our runs. One can 
prove 157 that Ei? n {t) > -T s d?/ ™(i), a negative number. In all the runs we 
considered we found E^? n > 0, always pulling up the negative ground-state 
value, but whether this is generally true, remains to be proven. 

This graph also demonstrates the importance of memory effects (and 
now our poor ground-stater is really running). In the top panel is a plot of 
the value of k e ff(t) described earlier. This parameter completely identifies 
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Fig. 5. Non-locality in time: the top panel shows k e ff(t), middle panel E c , bottom 
panel -Ec- The parameters in Eq.56 are k = 0.25, ui = 0.75 and e = 0.1. 



the density profile. The figure suggests that the correlation potential v c (t) 
is a highly non-local functional of the density, meaning that v c (t) depends 
not just on the density at and around time t, but rather on its entire history. 
The density profiles for a time range centered near t = 4.8 and centered 
near t = 28.9 are almost the same, yet the values of E c (t) near those times 
are hugely different. The density at times near t is not enough to specify E c : 
it depends on the entire history and this is what we mean by non-locality 
in time. Now E c (t) is intimately related to the correlation potential 92 : 



E c (t) = J d 3 rv c (rt)h(rt) 



(57) 



so that non-locality in E c directly implies non-locality in the correlation 
potential v c (t). Clearly, E c (t) will also be a highly non-local functional of 
the density and this is also shown in the figure. Any adiabatic approxima- 
tion has no memory and will fail to capture this effect. On the other hand, 
it may approximate the exchange potential well since, even for more than 
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two electrons, this is probably not a strongly non-local in time functional 
of the density 157 . 

We end this section by looking back to the linear response theory of 
the previous section. One can prove that although E^ vn vanishes in linear 
response, vt m , T* yn and V* vn do not 157 . 

2.8. Can we always find a Kohn-Sham potential for an 
appropriate initial state? 

The Runge-Gross theorem 1 proves that there is a 1-1 mapping between 
potentials and densities, much like in the ground-state case, but with one 
major difference: in the time-dependent case, the mapping is unique only 
for a specified initial state. There may be several different initial wave- 
functions that evolve with the same time-dependent density in different 
time-dependent potentials; the exchange-correlation potential for each of 
them will be different. All functional in use today completely ignore this 
initial-state dependence, partly because very little is known about it. 

To tackle the problem of how functional depend on the initial state we 
might begin by considering the simplest possible case: one electron. We ask, 
can we find two different initial states which evolve with the same density 
for all time in two different potentials? In fact, we cannot 159 ! There is 
no initial-state dependence for one electron. The proof is very simple: the 
two candidate wavefunctions must be related to each other by a phase, 
4>{rt) = e' a ( ri ' 4>{rt) in order to have the same density. Imposing equal 
h(rt) = — V • j(ri), gives V • [nVa] = 0. This can only be satisfied if a 
is a physically irrelevant constant. So in the one electron case, there is at 
most one wavefunction which can evolve with a given density: the evolving 
density uniquely specifies the potential. 

The situation is however quite different for two or more electrons. Given 
two initial wavefunctions that have the same density and first-time deriva- 
tive of the density 160 , and which are well-behaved in that their expectation 
values of the momentum-stress tensor are finite 159 , there exist two differ- 
ent potentials in which they each evolve with the same density for all time 
160 . This holds also in the case of two different inter-particle interactions 
and so is a statement about non-interacting w-representability. It was also 
shown 160 how to construct the difference in the potentials that keeps the 
two wavefunctions evolving with the same density. Essentially the proof 
follows from requiring the second derivative of the density to be the same 
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for each wavefunction, using the continuity equation to write this in terms 
of the current and then considering the Heisenberg equation of motion for 
the current 160 . This gives the first term in a Taylor series in time of the 
potential. The higher order terms may be obtained by considering higher 
order time-derivatives of the density; one ends up having to evaluate nested 
commutators of the momentum-stress tensor with the Hamiltonian. 

A simple example of initial-state dependence 159 occurs for two non- 
interacting electrons in one dimension. In Fig. 6 we give another example. 
The reference wavefunction is the lowest triplet state of the harmonic os- 
cillator, with orbitals shown in the figure. The orbitals of the other wave- 
function are also shown in the top panel; these are chosen so that the 
two wavefunctions have the same initial density and current. The reference 
wavefunction evolves with constant density in the harmonic oscillator po- 
tential shown below. The initial potential in which the second wavefunction 
evolves so that its density also remains constant is also shown. This alter- 
native potential will change in the next instant of time in order to keep 
the alternate orbitals evolving with the same constant density as the ref- 
erence ones. If we consider the density in the top panel to be that of an 
interacting two-electron system in some external potential, then the differ- 
ence between the two potentials in the lower panel is the difference in the 
exchange-correlation potential when the two different initial Kohn-Sham 
states of the top panel are chosen. 

2.9. What is memory? 

One of the first things a budding quantum mechanic learns is that knowl- 
edge of the many-electron wavefunction at any instant of time is enough to 
completely determine all properties of the system at that time. The Runge- 
Gross theorem 1 says that this is overkill: knowing just the time-evolving 
density and just the initial wavefunction is enough to know everything about 
the system. Like in the ground-state case, by trading in the 6A?"-variabled 
wavefunction for the 3-variabled density as its main player, and having to 
solve N uncoupled 3-variabled Schrodinger's equations instead of the daunt- 
ing coupled Schrodinger equation in 6N variables, TDDFT provides a more 
feasible framework for the dynamics. The subtleties of the electron-electron 
interplay are hidden in the functionals, which are in practice approximated. 

A feature of the time-dependent functionals is memory dependence, 
about which much has yet to be understood. Functionals in TDDFT are 
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Fig. 6. Initial-state dependence: the solid lines in the top panel are the two orbitals 
of one wavefunction, which happens to be stationary state of the harmonic oscillator 
shown below as a solid line. The dashed lines are the orbitals of an alternative initial 
wavefunction which evolve with the same density (thick solid line in top figure) in the 
potential which at time is shown as the dashed line in the bottom figure. 

haunted by the past: in general they depend on the density along its entire 
history, on the initial state of the interacting system and also on the choice 
of Kohn-Sham initial state. Examples of the history dependence and of the 
initial-state dependence have been given earlier in this chapter (Sees. 2.7 
and 2.8). 

How does this memory arise? In any wavefunction theory, no memory is 
needed: as our budding quantum mechanic would tell us, the wavefunction 
at time t contains all the information. The memory in TDDFT is a price we 
have to pay for the trading in of the complicated wavefunction for the much 
simpler density. In a sense, what we mean by memory is a consequence of 
the underlying (banished) wavefunction at time t not being a ground-state 
wavefunction. 

Recently it has been shown that the two sources of memory in TDDFT, 
history dependence and initial-state dependence, are inextricably inter- 
twined 154 . In fact, initial-state dependence can often be completely ab- 
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sorbed into a history-dependence. Imagine an interacting time-dependent 
calculation, run from £ = 0, generating a density n(rt). The exchange- 
correlation potential at time £ is determined by the density at all previous 
times as well as the initial interacting and Kohn-Sham states, and $o 
respectively. Now, if we knew what these wavefunctions evolved to at time 
£' < £, we may equally well think of £' as the initial time, and the inputs to 
v xc would be the density at times between £' and £ and the states , $r 
at time £'. This gives us an exact condition on the functional : 

w xc K (rt) = w X c[n;*o,*o](ri) for £ > £', (58) 

where 

n t >(Tt) = n(rt) for £ > t' (59) 

and rif (r£) is undefined for t < t' . Eq. (58) is a very difficult condition for 
approximate functionals to satisfy. Like other exact conditions, for example, 
the harmonic potential theorem 153 and the virial theorem 92 , it may be 
used as a feasibility test for approximate functionals. Any functional with 
history-dependence, and without initial-state dependence, very likely vio- 
lates this condition. On the other hand, this condition is trivially satisfied 
by any adiabatic approximation, which ignores both the dependence on the 
initial-state and on the history. 

An important consequence of Eq. (58) is that the initial-state depen- 
dence can often be completely expressed as a history effect along a "pseudo- 
prehistory" : once a system can be propagated backwards in time to some 
non-degenerate ground state, generating a density n(rt), then 

v XG [n; * , *o](rt) = v xc [n](rt) for t > 0. (60) 

Here, the density n(r£) = n(rt) for t > and is defined to be the density 
along a pseudo-prehistory which begins in some ground-state ^ gs ($ fls for 
the Kohn-Sham system) at some negative time t = t gs < and evolves 
under some many-electron Hamiltonian that carries us to the true ini- 
tial states ($o) at £ = 0. Initial-state dependence has vanished on the 
right-hand side of Eq. (60) since the systems "start" in the non-degenerate 
ground-states ^ gs and which, by the Hohenberg-Kohn theorem, are 
functionals of the ground-state density. Instead it has been absorbed into 
a pseudo-prehistory. 

When can we find a pseudo-prehistory and thus eliminate initial-state 
dependence? It can be shown that an arbitrary initial state cannot be 
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evolved back to some ground-state under a many-electron Hamiltonian 154 ; 
but there are plenty of initial states that can be. 

3. Conclusions: Dante's Inferno? 

When we begin applying TDDFT to quantum mechanical systems, the 
linear-response method looks very appealing, especially employing ALDA 
for /xc- We are like Dante at the edge of the woods. But once we study 
things a little more closely, we begin our descent through the various circles 
of Hell. So long as we stay within linear response, we are always looking at 
variations away from a ground state, and so we can consider ourselves in 
the outer circle. But once we begin to study fully time-dependent problems, 
then the initial-state dependence rears its ugly head (similar to that of 
Bertrand de Born) and we know we are in the inner circle. 

But, there is more than just the Inferno to Dante's classic. There follows 
the purgatorio and finally paradiso. It is necessary to first pass through Hell, 
to achieve the wisdom needed to reach the paradiso. For TDDFT, we are 
on our way. 
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